Imports
library(evd);
library(evdbayes);
library(coda);
library(ismev);
library(xts);
Loading the Data
load("../data/CAPE_Minder_Rychener_Malsot.RData")
load("../data/NINO34.RData")
load("../data/SRH_Minder_Rychener_Malsot.RData")
Generate PROD
prod = sqrt(cape)*srh
** Create Time Series Objects **
dates <- seq.Date(as.Date("1979-1-1"), as.Date("2015-12-31"), by=1)
feb29ix <- format(as.Date(dates), "%m-%d") == "02-29"
dates <- dates[!feb29ix]
prod_ts = xts(prod, order.by = rep(dates, each=8))
Beginning the Analysis
month_names = c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
get_monthly = function(x) {
output = list()
len = nrow(x)
# Get month of element
month = time(x)
month = gsub("....-", "", month)
month = gsub("-..", "", month)
monthlist = unique(month)
for (i in 1:12) {
output[[i]] = x[month == monthlist[i],]
}
names(output) = month_names
return(output)
}
monthly_max = get_monthly(apply.monthly(prod_ts, max))
r = 2
r_monthly_max = get_monthly(apply.monthly(prod_ts, function(x) order(x, decreasing=T)[1:r]))
# get the monthly maxima of prod
m1 = as.data.frame(apply.monthly(prod_ts, max))$V1;
# produce the gumbel qq plot
gumbel_qq = function(x, title="") qqplot(qgumbel(c(1:length(x))/(length(x)+1)),
x,
main = paste("Gumbell Q-Q Plot", title),
xlab = "Theoretical Quantiles",
ylab = "Sample Quantiles") ;
gumbel_qq(m1)

#qqplot(qgumbel(c(1:length(m1))/(length(m1)+1)),m1,main = "Gumbell Q-Q Plot",xlab = "Theoretical Quantiles", ylab = "Sample Quantiles") ;
The qq plot doesn’t fit very well, especially in the lower tail. This is likely due to seasonal dependence.
Fitting GEV to the entire data
# fit gevd with MLE and produce diagnostic plots
fitmax.MLE<-fgev(m1,method="Nelder-Mead")
par(mfrow=c(2,2))
fitmax.MLE
Call: fgev(x = m1, method = "Nelder-Mead")
Deviance: 9272.599
Estimates
loc scale shape
7.519e+03 7.013e+03 4.757e-03
Standard Errors
loc scale shape
421.90201 331.43749 0.05347
Optimization Information
Convergence: successful
Function Evaluations: 171
plot(fitmax.MLE)

Poor fit, probably because the distribution isn’t stationary. This is especially visible in the Probability plot, in which the confidence band is surpassed, indicating a poor fit.
# fit gevd with Bayesian Techniques
# use the previous outputs (rounded) as initial values (use a different shape)
init<-c(1.6e4,4e3,0.1)
# arbitrary priors
mat <- diag(c(10000,10000,100))
pn <- prior.norm(mean=c(0,0,0),cov=mat)
# proposal standard deviation (found by trial and error to get 20%<acceptance rate<40%)
psd<-c(500,0.03,0.02)
# draw 3k samples from markov chain
nit = 30000
thinning = 300
post<-posterior(nit, init=init,prior=pn,lh="gev",data=m1,psd=psd)
# diagnostic plots
MCMC<-mcmc(post[1:nit %% thinning == 0, ], thin=300)
plot(MCMC)

attr(mcmc(post),"ar")
mu sigma xi total
acc.rates 0.24 0.72 0.70 0.55
ext.rates 0.00 0.00 0.01 0.00
#MCMC_stab <- mcmc(post, thin=50, start=1000)
acf(mcmc(post[1:nit %% thinning == 0, ]))

We observe that there seem to be no substantial issues from the autocorrelation plots.
apply(mcmc(post[1:nit %% thinning == 0, ]),2,mean)
mu sigma xi
216.8681802 11212.0488049 -0.1898789
apply(mcmc(post[1:nit %% thinning == 0, ]),2,sd)
mu sigma xi
101.68289973 737.97773473 0.04013345
** Fit with MLE for months separately**
#monthly_fits = lapply(monthly_max,
# function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
monthly_fits = list()
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
print(i)
if (i %in% error_cases) {
monthly_fits[[i]] = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
std.err = FALSE)
}
else {
monthly_fits[[i]] = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead")
}
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
optimization may not have succeeded
[1] 10
[1] 11
[1] 12
names(monthly_fits) = names(monthly_max)
gumbel_qq(data.frame(monthly_max[[9]])[,1], "September")

gumbel_qq(data.frame(monthly_max[[12]])[,1], "December")

get_se = function(x, ix) {
if (is.null(x$std.err)) 0.01
else x$std.err[ix]
}
mle_loc = unlist(lapply(monthly_fits, function(x) x$estimate[1]))
mle_loc_se = unlist(lapply(monthly_fits, get_se, 1))
mle_scale = unlist(lapply(monthly_fits, function(x) x$estimate[2]))
mle_scale_se = unlist(lapply(monthly_fits, get_se, 2))
mle_shape = unlist(lapply(monthly_fits, function(x) x$estimate[3]))
mle_shape_se = unlist(lapply(monthly_fits, get_se, 3))
plot_w_err = function(x, y, se, title = NULL) {
max_ix = which.max(y)
min_ix = which.min(y)
plot(x, y,
ylim = c(y[min_ix] - se[min_ix], y[max_ix] + se[max_ix]),
main = title)
arrows(x,y-se,x,y+se, code=3, length=0.02, angle = 90)
}
plot_w_err(1:12, mle_loc, mle_loc_se, "Location Parameter as Estimated with Likelihood")

plot_w_err(1:12, mle_scale, mle_scale_se, "Scale Parameter as Estimated with Likelihood")

plot_w_err(1:12, mle_shape, mle_shape_se, "Shape Parameter as Estimated with Likelihood")

** Fit with Bayesian Methods for months separately**
bayes_fitter = function(x,
init = c(1e3, 1e3, 0.1), # Initial values
mat = diag(c(10000,10000,100)),
psd = c(500,0.1,0.1), # Proposed SDev
nit = 3000, # Nb Iterations
thinning = 50, # Thinning Factor
do_plots = FALSE, # Bool whether to show plot
seed = 42 # Seed
) {
set.seed(seed)
pn = prior.norm(mean=c(0,0,0),cov=mat)
post = posterior(nit, init=init, prior=pn, lh="gev", data=x, psd=psd)
if(do_plots) {
MCMC<-mcmc(post[1:nit %% thinning == 0, ])
plot(MCMC)
}
list(posterior = post,
acceptance_rate = attr(mcmc(post),"ar"))
}
monthly_bayes_fit = lapply(monthly_max, bayes_fitter, do_plots = TRUE,
psd = c(500,0.3,0.3), nit=30000, thinning = 300)












acceptance_rates = lapply(monthly_bayes_fit, function(x) x$acceptance_rate[1,])
print(acceptance_rates)
$jan
mu sigma xi total
0.16 0.56 0.58 0.44
$feb
mu sigma xi total
0.24 0.54 0.44 0.40
$mar
mu sigma xi total
0.24 0.32 0.20 0.25
$apr
mu sigma xi total
0.23 0.19 0.15 0.19
$may
mu sigma xi total
0.24 0.17 0.11 0.17
$jun
mu sigma xi total
0.24 0.18 0.12 0.18
$jul
mu sigma xi total
0.23 0.12 0.09 0.14
$aug
mu sigma xi total
0.24 0.22 0.14 0.20
$sep
mu sigma xi total
0.23 0.15 0.11 0.16
$oct
mu sigma xi total
0.24 0.38 0.23 0.28
$nov
mu sigma xi total
0.24 0.46 0.33 0.34
$dec
mu sigma xi total
0.19 0.58 0.53 0.43
rr # fit the r largest method #rl = rlarg.fit(data) data(venice) venfit <- rlarg.fit(venice[,-1])
$conv
[1] 0
$nllh
[1] 1139.09
$mle
[1] 120.5479027 12.7840265 -0.1129418
$se
[1] 1.36234055 0.54944881 0.01986948
PART 2
rr #plot the maximum as a function of ENSO plot(nino34,m1,main = of Monthly Maxima of PROD vs ENSO,xlab = , ylab = Maxima of PROD)

rr n_years = length(m1)/12 month_indexes = rep(c(1,2,3,4,5,6,7,8,9,10,11,12),n_years) plot(month_indexes,m1,main = Maxima of PROD,xlab = , ylab = Maxima of PROD)

# plot the chi plot for dependance analysis
dat.m1_month = cbind(m1,month_indexes);
dat.m1_nino = cbind(m1,nino34);
chiplot(dat.m1_month);
chiplot(dat.m1_nino);
PART 3 We will now analyse temporal clustering of extremes. For this, we will use the exiplot function from the evd library.
# get sample quantiles (used as plotting thresholds)
prod.quantiles = quantile(prod, c(0.8, 0.85, 0.9, 0.95,0.99,0.999))
exiplot(prod,prod.quantiles,main = "Analyiss of Extremal clustering",xlab = "Threshold", ylab = "Extremal Index")
We observe that the extremal index is ~0.2, we can therefore conclude that we have strong dependance of extremes, with the limiting mean cluster size being roughly 5. The clustering has no effect for estimaters based on the (monthly) maximum, but the r largest estimator is influenced by it.
LS0tCnRpdGxlOiAiVGh1bmRlcnN0cm9tIEFuYWx5c2lzIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCioqSW1wb3J0cyoqCmBgYHtyfQpsaWJyYXJ5KGV2ZCk7IApsaWJyYXJ5KGV2ZGJheWVzKTsKbGlicmFyeShjb2RhKTsKbGlicmFyeShpc21ldik7CmxpYnJhcnkoeHRzKTsKCmBgYAoKCioqTG9hZGluZyB0aGUgRGF0YSoqCgpgYGB7cn0KbG9hZCgiLi4vZGF0YS9DQVBFX01pbmRlcl9SeWNoZW5lcl9NYWxzb3QuUkRhdGEiKQpsb2FkKCIuLi9kYXRhL05JTk8zNC5SRGF0YSIpCmxvYWQoIi4uL2RhdGEvU1JIX01pbmRlcl9SeWNoZW5lcl9NYWxzb3QuUkRhdGEiKQpgYGAKCioqR2VuZXJhdGUgUFJPRCoqCmBgYHtyfQpwcm9kID0gc3FydChjYXBlKSpzcmgKYGBgCgoKCgoqKiBDcmVhdGUgVGltZSBTZXJpZXMgT2JqZWN0cyAqKgpgYGB7cn0KZGF0ZXMgPC0gc2VxLkRhdGUoYXMuRGF0ZSgiMTk3OS0xLTEiKSwgYXMuRGF0ZSgiMjAxNS0xMi0zMSIpLCBieT0xKQpmZWIyOWl4IDwtIGZvcm1hdChhcy5EYXRlKGRhdGVzKSwgIiVtLSVkIikgPT0gIjAyLTI5IgpkYXRlcyA8LSBkYXRlc1shZmViMjlpeF0KCnByb2RfdHMgPSB4dHMocHJvZCwgb3JkZXIuYnkgPSByZXAoZGF0ZXMsIGVhY2g9OCkpCmBgYAoKCgoqKkJlZ2lubmluZyB0aGUgQW5hbHlzaXMqKgpgYGB7cn0KbW9udGhfbmFtZXMgPSBjKCJqYW4iLCJmZWIiLCJtYXIiLCJhcHIiLCJtYXkiLCJqdW4iLCJqdWwiLCJhdWciLCJzZXAiLCJvY3QiLCJub3YiLCJkZWMiKQpnZXRfbW9udGhseSA9IGZ1bmN0aW9uKHgpIHsKICBvdXRwdXQgPSBsaXN0KCkKICBsZW4gPSBucm93KHgpCiAgCiAgIyBHZXQgbW9udGggb2YgZWxlbWVudAogIG1vbnRoID0gdGltZSh4KQogIG1vbnRoID0gZ3N1YigiLi4uLi0iLCAiIiwgbW9udGgpCiAgbW9udGggPSBnc3ViKCItLi4iLCAiIiwgbW9udGgpCiAgbW9udGhsaXN0ID0gdW5pcXVlKG1vbnRoKQogIGZvciAoaSBpbiAxOjEyKSB7CiAgICBvdXRwdXRbW2ldXSA9IHhbbW9udGggPT0gbW9udGhsaXN0W2ldLF0KICB9CiAgbmFtZXMob3V0cHV0KSA9IG1vbnRoX25hbWVzCiAgcmV0dXJuKG91dHB1dCkKfQoKbW9udGhseV9tYXggPSBnZXRfbW9udGhseShhcHBseS5tb250aGx5KHByb2RfdHMsIG1heCkpCnIgPSAyCnJfbW9udGhseV9tYXggPSBnZXRfbW9udGhseShhcHBseS5tb250aGx5KHByb2RfdHMsIGZ1bmN0aW9uKHgpIG9yZGVyKHgsIGRlY3JlYXNpbmc9VClbMTpyXSkpCmBgYAoKCmBgYHtyfQojIGdldCB0aGUgbW9udGhseSBtYXhpbWEgb2YgcHJvZAptMSA9IGFzLmRhdGEuZnJhbWUoYXBwbHkubW9udGhseShwcm9kX3RzLCBtYXgpKSRWMTsKIyBwcm9kdWNlIHRoZSBndW1iZWwgcXEgcGxvdApndW1iZWxfcXEgPSBmdW5jdGlvbih4LCB0aXRsZT0iIikgcXFwbG90KHFndW1iZWwoYygxOmxlbmd0aCh4KSkvKGxlbmd0aCh4KSsxKSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtYWluID0gcGFzdGUoIkd1bWJlbGwgUS1RIFBsb3QiLCB0aXRsZSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeGxhYiA9ICJUaGVvcmV0aWNhbCBRdWFudGlsZXMiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5bGFiID0gIlNhbXBsZSBRdWFudGlsZXMiKSA7CgpndW1iZWxfcXEobTEpCgojcXFwbG90KHFndW1iZWwoYygxOmxlbmd0aChtMSkpLyhsZW5ndGgobTEpKzEpKSxtMSxtYWluID0gIkd1bWJlbGwgUS1RIFBsb3QiLHhsYWIgPSAiVGhlb3JldGljYWwgUXVhbnRpbGVzIiwgeWxhYiA9ICJTYW1wbGUgUXVhbnRpbGVzIikgOwpgYGAKVGhlIHFxIHBsb3QgZG9lc24ndCBmaXQgdmVyeSB3ZWxsLCBlc3BlY2lhbGx5IGluIHRoZSBsb3dlciB0YWlsLiBUaGlzIGlzIGxpa2VseSBkdWUKdG8gc2Vhc29uYWwgZGVwZW5kZW5jZS4KCioqRml0dGluZyBHRVYgdG8gdGhlIGVudGlyZSBkYXRhKioKYGBge3J9CiMgZml0IGdldmQgd2l0aCBNTEUgYW5kIHByb2R1Y2UgZGlhZ25vc3RpYyBwbG90cwpmaXRtYXguTUxFPC1mZ2V2KG0xLG1ldGhvZD0iTmVsZGVyLU1lYWQiKQpwYXIobWZyb3c9YygyLDIpKQpmaXRtYXguTUxFCnBsb3QoZml0bWF4Lk1MRSkKYGBgClBvb3IgZml0LCBwcm9iYWJseSBiZWNhdXNlIHRoZSBkaXN0cmlidXRpb24gaXNuJ3Qgc3RhdGlvbmFyeS4gVGhpcyBpcyBlc3BlY2lhbGx5IAp2aXNpYmxlIGluIHRoZSBQcm9iYWJpbGl0eSBwbG90LCBpbiB3aGljaCB0aGUgY29uZmlkZW5jZSBiYW5kIGlzIHN1cnBhc3NlZCwgCmluZGljYXRpbmcgYSBwb29yIGZpdC4KCgpgYGB7cn0KIyBmaXQgZ2V2ZCB3aXRoIEJheWVzaWFuIFRlY2huaXF1ZXMKIyB1c2UgdGhlIHByZXZpb3VzIG91dHB1dHMgKHJvdW5kZWQpIGFzIGluaXRpYWwgdmFsdWVzICh1c2UgYSBkaWZmZXJlbnQgc2hhcGUpCmluaXQ8LWMoMS42ZTQsNGUzLDAuMSkKIyBhcmJpdHJhcnkgcHJpb3JzCm1hdCA8LSBkaWFnKGMoMTAwMDAsMTAwMDAsMTAwKSkgCnBuIDwtIHByaW9yLm5vcm0obWVhbj1jKDAsMCwwKSxjb3Y9bWF0KQojIHByb3Bvc2FsIHN0YW5kYXJkIGRldmlhdGlvbiAoZm91bmQgYnkgdHJpYWwgYW5kIGVycm9yIHRvIGdldCAyMCU8YWNjZXB0YW5jZSByYXRlPDQwJSkKcHNkPC1jKDUwMCwwLjAzLDAuMDIpCiMgZHJhdyAzayBzYW1wbGVzIGZyb20gbWFya292IGNoYWluCm5pdCA9IDMwMDAwCnRoaW5uaW5nID0gMzAwCnBvc3Q8LXBvc3RlcmlvcihuaXQsIGluaXQ9aW5pdCxwcmlvcj1wbixsaD0iZ2V2IixkYXRhPW0xLHBzZD1wc2QpCiMgZGlhZ25vc3RpYyBwbG90cwpNQ01DPC1tY21jKHBvc3RbMTpuaXQgJSUgdGhpbm5pbmcgPT0gMCwgXSwgdGhpbj0zMDApIApwbG90KE1DTUMpIAphdHRyKG1jbWMocG9zdCksImFyIikKCmBgYAoKCmBgYHtyfQojTUNNQ19zdGFiIDwtIG1jbWMocG9zdCwgdGhpbj01MCwgc3RhcnQ9MTAwMCkKYWNmKG1jbWMocG9zdFsxOm5pdCAlJSB0aGlubmluZyA9PSAwLCBdKSkKYGBgCldlIG9ic2VydmUgdGhhdCB0aGVyZSBzZWVtIHRvIGJlIG5vIHN1YnN0YW50aWFsIGlzc3VlcyBmcm9tIHRoZSBhdXRvY29ycmVsYXRpb24gCnBsb3RzLiAKCmBgYHtyfQphcHBseShtY21jKHBvc3RbMTpuaXQgJSUgdGhpbm5pbmcgPT0gMCwgXSksMixtZWFuKQphcHBseShtY21jKHBvc3RbMTpuaXQgJSUgdGhpbm5pbmcgPT0gMCwgXSksMixzZCkKCmBgYAoKKiogRml0IHdpdGggTUxFIGZvciBtb250aHMgc2VwYXJhdGVseSoqCmBgYHtyfQojbW9udGhseV9maXRzID0gbGFwcGx5KG1vbnRobHlfbWF4LCAKIyAgICAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSBmZ2V2KGRhdGEuZnJhbWUoeClbLDFdLCBtZXRob2Q9Ik5lbGRlci1NZWFkIikpCm1vbnRobHlfZml0cyA9IGxpc3QoKQplcnJvcl9jYXNlcyA9IGMoOSwgMTIpCmZvciAoaSBpbiAxOmxlbmd0aChtb250aGx5X21heCkpIHsKICBwcmludChpKQogIAogIGlmIChpICVpbiUgZXJyb3JfY2FzZXMpIHsKICAgIG1vbnRobHlfZml0c1tbaV1dID0gZmdldihhcy5kYXRhLmZyYW1lKG1vbnRobHlfbWF4W1tpXV0pJFYxLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtZXRob2QgPSAiTmVsZGVyLU1lYWQiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0ZC5lcnIgPSBGQUxTRSkKICB9CiAgZWxzZSB7CiAgICBtb250aGx5X2ZpdHNbW2ldXSA9IGZnZXYoYXMuZGF0YS5mcmFtZShtb250aGx5X21heFtbaV1dKSRWMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gIk5lbGRlci1NZWFkIikKICB9Cn0KbmFtZXMobW9udGhseV9maXRzKSA9IG5hbWVzKG1vbnRobHlfbWF4KQoKYGBgCgpgYGB7cn0KZ3VtYmVsX3FxKGRhdGEuZnJhbWUobW9udGhseV9tYXhbWzldXSlbLDFdLCAiU2VwdGVtYmVyIikKZ3VtYmVsX3FxKGRhdGEuZnJhbWUobW9udGhseV9tYXhbWzEyXV0pWywxXSwgIkRlY2VtYmVyIikKYGBgCmBgYHtyfQpnZXRfc2UgPSBmdW5jdGlvbih4LCBpeCkgewogIGlmIChpcy5udWxsKHgkc3RkLmVycikpIDAKICBlbHNlIHgkc3RkLmVycltpeF0KfQptbGVfbG9jID0gdW5saXN0KGxhcHBseShtb250aGx5X2ZpdHMsIGZ1bmN0aW9uKHgpIHgkZXN0aW1hdGVbMV0pKQptbGVfbG9jX3NlID0gdW5saXN0KGxhcHBseShtb250aGx5X2ZpdHMsIGdldF9zZSwgMSkpCm1sZV9zY2FsZSA9IHVubGlzdChsYXBwbHkobW9udGhseV9maXRzLCBmdW5jdGlvbih4KSB4JGVzdGltYXRlWzJdKSkKbWxlX3NjYWxlX3NlID0gdW5saXN0KGxhcHBseShtb250aGx5X2ZpdHMsIGdldF9zZSwgMikpCm1sZV9zaGFwZSA9IHVubGlzdChsYXBwbHkobW9udGhseV9maXRzLCBmdW5jdGlvbih4KSB4JGVzdGltYXRlWzNdKSkKbWxlX3NoYXBlX3NlID0gdW5saXN0KGxhcHBseShtb250aGx5X2ZpdHMsIGdldF9zZSwgMykpCmBgYAoKYGBge3J9CnBsb3Rfd19lcnIgPSBmdW5jdGlvbih4LCB5LCBzZSwgdGl0bGUgPSBOVUxMKSB7CiAgbWF4X2l4ID0gd2hpY2gubWF4KHkpCiAgbWluX2l4ID0gd2hpY2gubWluKHkpCiAgcGxvdCh4LCB5LAogICAgICAgeWxpbSA9IGMoeVttaW5faXhdIC0gc2VbbWluX2l4XSwgeVttYXhfaXhdICsgc2VbbWF4X2l4XSksCiAgICAgICBtYWluID0gdGl0bGUpCiAgYXJyb3dzKHgseS1zZSx4LHkrc2UsIGNvZGU9MywgbGVuZ3RoPTAuMDIsIGFuZ2xlID0gOTApCn0KcGxvdF93X2VycigxOjEyLCBtbGVfbG9jLCBtbGVfbG9jX3NlLCAiTG9jYXRpb24gUGFyYW1ldGVyIGFzIEVzdGltYXRlZCB3aXRoIExpa2VsaWhvb2QiKQpwbG90X3dfZXJyKDE6MTIsIG1sZV9zY2FsZSwgbWxlX3NjYWxlX3NlLCAiU2NhbGUgUGFyYW1ldGVyIGFzIEVzdGltYXRlZCB3aXRoIExpa2VsaWhvb2QiKQpwbG90X3dfZXJyKDE6MTIsIG1sZV9zaGFwZSwgbWxlX3NoYXBlX3NlLCAiU2hhcGUgUGFyYW1ldGVyIGFzIEVzdGltYXRlZCB3aXRoIExpa2VsaWhvb2QiKQoKYGBgCgoqKiBGaXQgd2l0aCBCYXllc2lhbiBNZXRob2RzIGZvciBtb250aHMgc2VwYXJhdGVseSoqCmBgYHtyfQoKYmF5ZXNfZml0dGVyID0gZnVuY3Rpb24oeCwgCiAgICAgICAgICAgICAgICAgICAgICAgIGluaXQgPSBjKDFlMywgMWUzLCAwLjEpLCAjIEluaXRpYWwgdmFsdWVzCiAgICAgICAgICAgICAgICAgICAgICAgIG1hdCA9IGRpYWcoYygxMDAwMCwxMDAwMCwxMDApKSwKICAgICAgICAgICAgICAgICAgICAgICAgcHNkID0gYyg1MDAsMC4xLDAuMSksICMgUHJvcG9zZWQgU0RldgogICAgICAgICAgICAgICAgICAgICAgICBuaXQgPSAzMDAwLCAjIE5iIEl0ZXJhdGlvbnMKICAgICAgICAgICAgICAgICAgICAgICAgdGhpbm5pbmcgPSA1MCwgIyBUaGlubmluZyBGYWN0b3IKICAgICAgICAgICAgICAgICAgICAgICAgZG9fcGxvdHMgPSBGQUxTRSwgIyBCb29sIHdoZXRoZXIgdG8gc2hvdyBwbG90CiAgICAgICAgICAgICAgICAgICAgICAgIHNlZWQgPSA0MiAjIFNlZWQgCiAgICAgICAgICAgICAgICAgICAgICAgICkgewogIHNldC5zZWVkKHNlZWQpICAgICAgICAgICAgICAgIAogIHBuID0gcHJpb3Iubm9ybShtZWFuPWMoMCwwLDApLGNvdj1tYXQpCiAgcG9zdCA9IHBvc3RlcmlvcihuaXQsIGluaXQ9aW5pdCwgcHJpb3I9cG4sIGxoPSJnZXYiLCBkYXRhPXgsIHBzZD1wc2QpCiAgCiAgaWYoZG9fcGxvdHMpIHsKICAgIE1DTUM8LW1jbWMocG9zdFsxOm5pdCAlJSB0aGlubmluZyA9PSAwLCBdKSAKICAgIHBsb3QoTUNNQykgCiAgfQogIGxpc3QocG9zdGVyaW9yID0gcG9zdCwgCiAgICAgICBhY2NlcHRhbmNlX3JhdGUgPSBhdHRyKG1jbWMocG9zdCksImFyIikpCn0KCm1vbnRobHlfYmF5ZXNfZml0ID0gbGFwcGx5KG1vbnRobHlfbWF4LCBiYXllc19maXR0ZXIsIGRvX3Bsb3RzID0gVFJVRSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgcHNkID0gYyg1MDAsMC4zLDAuMyksIG5pdD0zMDAwMCwgdGhpbm5pbmcgPSAzMDApCmFjY2VwdGFuY2VfcmF0ZXMgPSBsYXBwbHkobW9udGhseV9iYXllc19maXQsIGZ1bmN0aW9uKHgpIHgkYWNjZXB0YW5jZV9yYXRlWzEsXSkKcHJpbnQoYWNjZXB0YW5jZV9yYXRlcykKCmBgYAoKCgpgYGB7cn0KIyBmaXQgdGhlIHIgbGFyZ2VzdCBtZXRob2QKI3JsID0gcmxhcmcuZml0KGRhdGEpCmRhdGEodmVuaWNlKQp2ZW5maXQgPC0gcmxhcmcuZml0KHZlbmljZVssLTFdKQpgYGAKCgoqKlBBUlQgMioqCmBgYHtyfQojcGxvdCB0aGUgbWF4aW11bSBhcyBhIGZ1bmN0aW9uIG9mIEVOU08KcGxvdChuaW5vMzQsbTEsbWFpbiA9ICJDb21wYXJpc29uIG9mIE1vbnRobHkgTWF4aW1hIG9mIFBST0QgdnMgRU5TTyIseGxhYiA9ICJFTlNPIiwgeWxhYiA9ICJNb250aGx5IE1heGltYSBvZiBQUk9EIikKbl95ZWFycyA9IGxlbmd0aChtMSkvMTIKbW9udGhfaW5kZXhlcyA9IHJlcChjKDEsMiwzLDQsNSw2LDcsOCw5LDEwLDExLDEyKSxuX3llYXJzKQpwbG90KG1vbnRoX2luZGV4ZXMsbTEsbWFpbiA9ICJNb250aGx5IE1heGltYSBvZiBQUk9EIix4bGFiID0gIk1vbnRoIiwgeWxhYiA9ICJNb250aGx5IE1heGltYSBvZiBQUk9EIikKYGBgCmBgYHtyfQojIHBsb3QgdGhlIGNoaSBwbG90IGZvciBkZXBlbmRhbmNlIGFuYWx5c2lzCmRhdC5tMV9tb250aCA9IGNiaW5kKG0xLG1vbnRoX2luZGV4ZXMpOwpkYXQubTFfbmlubyA9IGNiaW5kKG0xLG5pbm8zNCk7CmNoaXBsb3QoZGF0Lm0xX21vbnRoKTsKY2hpcGxvdChkYXQubTFfbmlubyk7CmBgYAoKCgoqKlBBUlQgMyoqCldlIHdpbGwgbm93IGFuYWx5c2UgdGVtcG9yYWwgY2x1c3RlcmluZyBvZiBleHRyZW1lcy4gRm9yIHRoaXMsIHdlIHdpbGwgdXNlIHRoZSBleGlwbG90IGZ1bmN0aW9uIGZyb20gdGhlIGV2ZCBsaWJyYXJ5LgoKYGBge3J9CiMgZ2V0IHNhbXBsZSBxdWFudGlsZXMgKHVzZWQgYXMgcGxvdHRpbmcgdGhyZXNob2xkcykKcHJvZC5xdWFudGlsZXMgPSBxdWFudGlsZShwcm9kLCBjKDAuOCwgMC44NSwgMC45LCAwLjk1LDAuOTksMC45OTkpKQpleGlwbG90KHByb2QscHJvZC5xdWFudGlsZXMsbWFpbiA9ICJBbmFseWlzcyBvZiBFeHRyZW1hbCBjbHVzdGVyaW5nIix4bGFiID0gIlRocmVzaG9sZCIsIHlsYWIgPSAiRXh0cmVtYWwgSW5kZXgiKQpgYGAKCldlIG9ic2VydmUgdGhhdCB0aGUgZXh0cmVtYWwgaW5kZXggaXMgfjAuMiwgd2UgY2FuIHRoZXJlZm9yZSBjb25jbHVkZSB0aGF0IHdlIGhhdmUgc3Ryb25nIGRlcGVuZGFuY2Ugb2YgZXh0cmVtZXMsIHdpdGggdGhlIGxpbWl0aW5nIG1lYW4gY2x1c3RlciBzaXplIGJlaW5nIHJvdWdobHkgNS4gVGhlIGNsdXN0ZXJpbmcgaGFzIG5vIGVmZmVjdCBmb3IgZXN0aW1hdGVycyBiYXNlZCBvbiB0aGUgKG1vbnRobHkpIG1heGltdW0sIGJ1dCB0aGUgciBsYXJnZXN0IGVzdGltYXRvciBpcyBpbmZsdWVuY2VkIGJ5IGl0Lgo=